Muon Lifetime

Santiago R. , Birge Súkrú "Ken"

In [227]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as cs
import pandas as pd
Data_Path = "Data/"
In [228]:
100*1e3/(2*1e-6)/cs.c
Out[228]:
166.78204759907604
In [229]:
v = cs.c/10
v*2*1e-6
Out[229]:
59.9584916

Energy Losses

In [230]:
Beta = 0.997
Gamma = 1/np.sqrt(1-Beta**2)
z = 1 #muon charge
K = -4*np.pi/(cs.m_e*cs.c**2)
K_2 = (cs.e**2/(4*np.pi*cs.epsilon_0))**2
n_al = cs.N_A*13*2.7*1e6/(26.981539)
n_lead = cs.N_A*82*11.34*1e6/(207.2)
I_Al = 10*cs.e*13
I_Lead = 10*cs.e*82
ln_Al = 2*cs.m_e*cs.c**2*Beta**2/(I_Al*(1-Beta**2))
ln_Lead = 2*cs.m_e*cs.c**2*Beta**2/(I_Lead*(1-Beta**2))
dEdx_Al = K*n_al*z**2/Beta**2*K_2*(np.log(ln_Al)-Beta**2) *1e-6*1e-2/cs.e
dEdx_Lead = K*n_lead*z**2/Beta**2*K_2*(np.log(ln_Lead)-Beta**2) *1e-6*1e-2/cs.e
print(dEdx_Al, "MeV/cm")
print(dEdx_Lead, "MeV/cm")

-5.259466094178444 MeV/cm
-15.590779410671892 MeV/cm
In [231]:
dEdx_Al*6
Out[231]:
-31.556796565070666
In [232]:
dEdx_Lead*10
Out[232]:
-155.90779410671894
In [233]:
1.6726219*1e-27*cs.c*5/cs.e
Out[233]:
15.648693785208149

On-Site Measurements and calculations

Runtime behavior of signals

Event Rates - Photomultipliers, Light Room

In [234]:
#PMT1
dNdt = np.array([713,673,670,673,667,670,681,742,639,651])
dt = 1 #second
Impedanz = 50 #ohm, Quad Scaler
In [235]:
#PMT2
dNdt = np.array([1390,1282,1470,1370,1301,1250,1201,1238,1244,1547])
dt = 1 #second
In [236]:
#PMT3
dNdt = np.array([44712,42164,39170,37567,39240,38539,36668,36473,39693,36576])
dt = 1 #second
In [237]:
#PMT4
dNdt = np.array([829,885,838,879,836,812,894,836,899,865])
dt = 1 #second
In [238]:
#PMT5
dNdt = np.array([23246,22930,23513,19091,19692,20946,20712,22777,20511,20160])
dt = 1 #second

Event Rates - Photomultipliers, Dark Room

In [239]:
#PMT1
dNdt = np.array([508,456,422,488,507,481,532,534,580,492])
dt = 1 #second
Impedanz = 50 #ohm, Quad Scaler
In [240]:
#PMT2
dNdt = np.array([1335,1360,1202,1278,1332,1163,1271,1338,1308,1317])
dt = 1 #second
In [241]:
#PMT3
dNdt = np.array([9502,9084,8716,9737,6833,6973,6758,6815,6656,6544,5417,9794])
dt = 1 #second
In [242]:
#PMT4
dNdt = np.array([726,901,834,924,982,951,858,854,851,945,767,800])
dt = 1 #second
In [243]:
#PMT5
dNdt = np.array([5414,6171,5245,7166,5224,5454,7847,7833,8470,9521,6790,6401])
dt = 1 #second

Event Rates - Scintillators

In [244]:
#SC1 - same as PMT 1
dNdt = np.array([508,456,422,488,507,481,532,534,580,492])
dt = 1 #second
Impedanz = 50 #ohm, Quad Scaler
In [245]:
#SC2 - PMT 2 and 4
dNdt = np.array([153,147,144,144,135,119,165,165,130,163])
dt = 1 #second
In [246]:
#SC3 - PMT 3 and 5
dNdt = np.array([156,170,171,157,149,162,129,172,161,173,164,154])
dt = 1 #second
In [247]:
#SC1 and SC2
dNdt = np.array([19,14,14,11,19,11,8,11,11,16,12,27])
dt = 1 #second
In [248]:
#SC1, SC2 and SC3
dNdt = np.array([2,5,8,4,3,4,10,4,8,6,7,2])
dt = 1 #second
In [249]:
#SC1 and SC2 with SC3 Veto
dNdt = np.array([11,6,8,11,15,13,8,12,13,13])
dt = 1 #second
In [250]:
#SC2 and SC3
dNdt = np.array([16,20,15,20,22,17,14,14,12,25])
dt = 1 #second
In [251]:
#SC2 or SC3
dNdt = np.array([348,302,304,261,306,288,301,287,291,269])
dt = 1 #second

Data Reading

In [252]:
Data = pd.read_csv(Data_Path+"MainRun1.Spe")
In [253]:
Data
Out[253]:
$SPEC_ID:
0 No sample description was entered.
1 $SPEC_REM:
2 DET# 1
3 DETDESC# FP03 Model 928 SN 19353133
4 AP# Maestro Version 7.01
... ...
16404 3
16405 -5.212607E+001 2.996516E+000 0.000000E+000 ns
16406 $SHAPE_CAL:
16407 3
16408 0.000000E+000 0.000000E+000 0.000000E+000

16409 rows × 1 columns

In [254]:
x = int(len(Data))-10
A = np.array(Data['$SPEC_ID:'][11:16000], dtype=int)
In [255]:
A
Out[255]:
array([0, 0, 0, ..., 0, 0, 0])
In [256]:
channels = np.arange(1,len(A)+1,1)
plt.plot(channels,A)
Out[256]:
[<matplotlib.lines.Line2D at 0x7f8ea9b010d0>]

Signal Delays

In [257]:
#Quad Coincidence Passthrough SC1 to Coincidence Unit Start Signal Logic
dt = 18 #ns

#Quad Coincidence SC2 or SC3 to Coincidence Unit Stop Signal Generator/Logic
dt = 15.6 #ns

Muon Lifetime [Dr. Marzieh Bahmani]

In [258]:
%matplotlib inline
import importlib
if not importlib.util.find_spec('iminuit'):
    !conda install -c conda-forge iminuit
import iminuit
%load_ext autoreload
%autoreload 1
%run muonutils

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

<Figure size 432x288 with 0 Axes>

Time Calibration Parser

In [259]:
import numpy as np
import matplotlib.pyplot as plt

with open("Data/Calibration/2mys.Spe") as fp:
    Lines = fp.readlines()
    print(len(Lines))
    count = 0
    for line in Lines:
        count += 1
        if "0 16383" in line:
            start = count
            break
    print("done")
    print(start)
    end = start+16383
    print(end)
    
x = np.arange(0,16383)

16410
done
12
16395
In [260]:
with open("Data/Calibration/80ns.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

120
Out[260]:
[<matplotlib.lines.Line2D at 0x7f8ea970c150>]
In [261]:
with open("Data/Calibration/2mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

682
Out[261]:
[<matplotlib.lines.Line2D at 0x7f8ea908bc10>]
In [262]:
with open("Data/Calibration/4_04.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

1367
Out[262]:
[<matplotlib.lines.Line2D at 0x7f8ea9df0350>]
In [263]:
with open("Data/Calibration/6_02mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

2031
Out[263]:
[<matplotlib.lines.Line2D at 0x7f8ea8e5b550>]
In [264]:
with open("Data/Calibration/8_04mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

2716
Out[264]:
[<matplotlib.lines.Line2D at 0x7f8ea95b4a50>]
In [265]:
with open("Data/Calibration/10mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

3387
Out[265]:
[<matplotlib.lines.Line2D at 0x7f8ea95e43d0>]
In [266]:
with open("Data/Calibration/12mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

4112
Out[266]:
[<matplotlib.lines.Line2D at 0x7f8ea8ea1a10>]
In [267]:
with open("Data/Calibration/12_5mys.Spe") as fp:
    Lines = fp.readlines()
    data = [int(i) for i in Lines[start:end]]
    dta = dta + np.array(data)

for i,val in enumerate(data):
    if val != 0:
        print(i+1)
        break
plt.plot(x,data,)

4242
Out[267]:
[<matplotlib.lines.Line2D at 0x7f8ea97c4b10>]
In [268]:
lines = ['       ' + str(i) + '\n' for i in dta]
with open("Data/Calibration/Calibration Parsed.Spe", "w") as fp:
    Lines[start:end] = lines
    fp.writelines(Lines)

plt.plot(x,dta)    
del dta

Background Signal and Delay

In [269]:
def draw_lifetime_spectrum(spectrum, fitresult=None):
    fig, ax = plt.subplots(figsize=(7., 5.))
    ax.set_ylabel('Entries')
    if spectrum.axislabel is not None:
        ax.set_xlabel(spectrum.axislabel)
    ax.margins(x=0)
    ax.hist(x=spectrum.centers, bins=spectrum.edges, weights=spectrum.values, histtype='step')
    if fitresult is not None:
        x = np.linspace(spectrum.edges[0], spectrum.edges[-1], 501)
        ax.plot(x, fitresult.curve(*fitresult.values)(x), zorder=0)
In [270]:
draw_lifetime_spectrum(load_spectrum("Data/Background.Spe"))
plt.savefig('background_spectrum.pdf')

Calibration

In [271]:
def draw_calibration_spectrum(spectrum, points):
    fig, ax = plt.subplots(nrows=2, sharex=True, gridspec_kw=dict(hspace=0), figsize=(7., 7.))
    ax[-1].set_xlabel('Channel number')
    ax[0].set_ylabel('Entries')
    ax[0].margins(x=0, y=0.25)
    ax[0].hist(bins=spectrum.edges, x=spectrum.centers, weights=spectrum.values, histtype='step')
    for i in range(len(points)):
        chan_no = points.at[i,'chan']
        text = '%.1f' % points.at[i,'time']
        color = 'gray'
        ax[0].annotate(text, xy=(chan_no, spectrum.values[chan_no]), xytext=(0., 15.), va='baseline', ha='center', fontsize='smaller', textcoords='offset points', color=color, arrowprops=dict(arrowstyle='->', shrinkA=0.1, facecolor=color, edgecolor=color))
    ax[1].set_ylabel('Time [μs]')
    ax[1].errorbar(points.chan, points.time, xerr=points.chan_err, yerr=points.time_err, fmt='none', edgecolor='k')

def draw_calibration_fit(points, fitresult):
    fig, ax = plt.subplots(figsize=(5., 5.))
    ax.set_xlabel('Channel number')
    ax.set_ylabel('Time [μs]')
    ax.errorbar(points.chan, points.time, xerr=points.chan_err, yerr=points.time_err, fmt='none', edgecolor='k', zorder=1)
    if fitresult is not None:
        x = np.linspace(0., 4296., 4296, False)
        ax.plot(x, fitresult.curve(*fitresult.values)(x), zorder=0)
In [272]:
calibration_points = pd.DataFrame(columns=['chan', 'time', 'chan_err', 'time_err'], data=[
# Format:
#       (Channel number, time, error of channel number, error of time)
# Beispiel:
       (120,   0.08,  5,  0.02),
       (682,  2.,  5,  0.2),
       (1367,  4.04,  5,  0.2),
       (2031,  6.02,  5,  0.2),
       (2716,  8.04,  5,  0.2),
       (3387,  10.,  5,  0.2),
       (4112,  12.,  5,  0.2),
       (4242,  12.5,  5,  0.2),
    ])
draw_calibration_spectrum(load_spectrum('Data/Calibration/Calibration Parsed.Spe'), calibration_points)
plt.savefig('calibration_spectrum.pdf')
In [273]:
def fit_calibration(points):
    def line(m, b):
        '''Return a line with slope m and intercept b.'''
        return lambda x: m * x + b
    return fit(method=weighted_least_squares,
               curve=line,
               data=dict(
                   x=points.chan.values,
                   y=points.time.values,
                   ey=points.time_err.values
               ),
               params=dict(
                   m=0.01, error_m=0.01,
                   b=0., error_b=0.01
               ))

calibration = fit_calibration(calibration_points)
In [274]:
calibration
Out[274]:
FCN = 3.643 Ncalls = 37 (56 total)
EDM = 1.23E-22 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 m 3.028E-3 0.027E-3
1 b -0.279 0.021
m b
m 1.00 -0.38
b -0.38 1.00
In [ ]:
In [275]:
draw_calibration_fit(calibration_points, calibration)
plt.savefig('calibration_fit.pdf')
In [276]:
calibration.values
Out[276]:
array([ 0.00302791, -0.27909678])
In [277]:
calibration.fval
Out[277]:
3.642617317060767
In [278]:
calibration.curve(*calibration.values)(np.array([500, 1000, 1500]))
Out[278]:
array([1.23485933, 2.74881545, 4.26277156])
In [279]:
3.643/5
Out[279]:
0.7285999999999999

Fit

In [280]:
raw_spectrum = load_spectrum('Data/MainRunFinal.Spe')
draw_lifetime_spectrum(raw_spectrum)
plt.savefig('raw_spectrum.pdf')
In [295]:
spectrum = raw_spectrum
spectrum = spectrum.rebin(start=1400, stop=3500, ngroup=10)
spectrum = spectrum.replace_axis(calibration.curve(*calibration.values)(spectrum.edges), 'Time [μs]')
draw_lifetime_spectrum(spectrum)
plt.savefig('rebinned_spectrum.pdf')
In [296]:
def fit_spectrum(spectrum, **params):
    # Define curve to be fitted:
    from numpy import exp
    def curve(tau_mu, n, c):
        f = 1.27
        Q = 0.993
        A = 0.7054
        tau_eff = 1/(Q/tau_mu+A)
        return lambda t: n / tau_mu * exp(-t/tau_mu) + n / (f*tau_eff) * exp(-t/tau_eff) + c
    # Default parameter configuration:
    params.setdefault('tau_mu', 1.)
    params.setdefault('error_tau_mu', params['tau_mu'])
    params.setdefault('n', 1000.)
    params.setdefault('error_n', params['n'])
    params.setdefault('c', spectrum.values[-4:].mean())
    params.setdefault('error_c', params['c'])
    # Run fit at the bin centers:
    return fit(method=chi_squared,
               curve=curve,
               data=spectrum,
               params=params)
In [297]:
fitresult = fit_spectrum(spectrum)
display(fitresult) 
draw_lifetime_spectrum(spectrum, fitresult)
plt.savefig('fit_spectrum.pdf')
FCN = 161.5 Ncalls = 166 (195 total)
EDM = 1.49E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 tau_mu 2.17 0.18
1 n 940 80
2 c 15.7 1.5
tau_mu n c
tau_mu 1.00 -0.95 -0.93
n -0.95 1.00 0.79
c -0.93 0.79 1.00

Systematic Errors

In [284]:
spectrum = raw_spectrum
spectrum = spectrum.rebin(start=1500, stop=3500, ngroup=10)
modified_calibration = calibration.values- calibration.values*(0.013/0.3)
spectrum = spectrum.replace_axis(calibration.curve(*modified_calibration)(spectrum.edges), 'Time [μs]')
draw_lifetime_spectrum(spectrum)
print(modified_calibration)

[ 0.0028967  -0.26700259]
In [285]:
fitresult = fit_spectrum(spectrum)
display(fitresult) 
draw_lifetime_spectrum(spectrum, fitresult)
FCN = 156.7 Ncalls = 170 (199 total)
EDM = 2.98E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 tau_mu 2.08 0.20
1 n 890 90
2 c 15.7 1.7
tau_mu n c
tau_mu 1.00 -0.96 -0.93
n -0.96 1.00 0.82
c -0.93 0.82 1.00
In [286]:
spectrum = raw_spectrum
spectrum = spectrum.rebin(start=1500, stop=3500, ngroup=10)
modified_calibration = calibration.values+ calibration.values*(0.013/0.3)
spectrum = spectrum.replace_axis(calibration.curve(*modified_calibration)(spectrum.edges), 'Time [μs]')
draw_lifetime_spectrum(spectrum)
print(modified_calibration)

[ 0.00315912 -0.29119098]
In [287]:
fitresult = fit_spectrum(spectrum)
display(fitresult) 
draw_lifetime_spectrum(spectrum, fitresult)
FCN = 156.7 Ncalls = 136 (165 total)
EDM = 1.38E-05 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 tau_mu 2.23 0.21
1 n 1000 110
2 c 15.8 1.6
tau_mu n c
tau_mu 1.00 -0.96 -0.93
n -0.96 1.00 0.82
c -0.93 0.82 1.00
In [288]:
#Different Binning
spectrum = raw_spectrum
spectrum = spectrum.rebin(start=1500, stop=8000, ngroup=20)
modified_calibration = calibration.values
spectrum = spectrum.replace_axis(calibration.curve(*modified_calibration)(spectrum.edges), 'Time [μs]')
draw_lifetime_spectrum(spectrum)
print(modified_calibration)

[ 0.00302791 -0.27909678]
In [289]:
fitresult = fit_spectrum(spectrum)
display(fitresult) 
draw_lifetime_spectrum(spectrum, fitresult)
FCN = 285.5 Ncalls = 139 (168 total)
EDM = 1.32E-06 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 tau_mu 2.74 0.09
1 n 1570 60
2 c 22.6 0.4
tau_mu n c
tau_mu 1.00 -0.90 -0.64
n -0.90 1.00 0.43
c -0.64 0.43 1.00
In [290]:
#Different Binning
spectrum = raw_spectrum
spectrum = spectrum.rebin(start=1500, stop=16000, ngroup=29)
modified_calibration = calibration.values
spectrum = spectrum.replace_axis(calibration.curve(*modified_calibration)(spectrum.edges), 'Time [μs]')
draw_lifetime_spectrum(spectrum)
print(modified_calibration)

[ 0.00302791 -0.27909678]
In [291]:
fitresult = fit_spectrum(spectrum)
display(fitresult) 
draw_lifetime_spectrum(spectrum, fitresult)
FCN = 487 Ncalls = 132 (159 total)
EDM = 1.47E-07 (Goal: 1E-05) up = 1.0
Valid Min. Valid Param. Above EDM Reached call limit
True True False False
Hesse failed Has cov. Accurate Pos. def. Forced
False True True True False
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 tau_mu 3.05 0.09
1 n 2100 70
2 c 29.77 0.30
tau_mu n c
tau_mu 1.00 -0.88 -0.44
n -0.88 1.00 0.27
c -0.44 0.27 1.00
In [ ]: